// Magic Software, Inc.
// http://www.magic-software.com
// Copyright (c) 2000, All Rights Reserved
//
// Source code from Magic Software is supplied under the terms of a license
// agreement and may not be copied or disclosed except in accordance with the
// terms of that agreement.  The various license agreements may be found at
// the Magic Software web site.  This file is subject to the license
//
// FREE SOURCE CODE
// http://www.magic-software.com/License/free.pdf

#ifndef MGCLINEARSYSTEM_H
#define MGCLINEARSYSTEM_H

#include "MgcMath.h"


class MgcLinearSystem
{
public:
    MgcLinearSystem () { /**/ }

    // 2x2 and 3x3 systems (avoids overhead of Gaussian elimination)
    static MgcReal& Tolerance ();
    static bool Solve2 (const MgcReal aafA[2][2], const MgcReal afB[2],
        MgcReal afX[2]);
    static bool Solve3 (const MgcReal aafA[3][3], const MgcReal afB[3],
        MgcReal afX[3]);

    // convenience for allocation, memory is zeroed out
    static MgcReal** NewMatrix (int iSize);
    static void DeleteMatrix (int iSize, MgcReal** aafA);
    static MgcReal* NewVector (int iSize);
    static void DeleteVector (int iSize, MgcReal* afB);

    static bool Inverse (int iSize, MgcReal** aafA);
    // Input:
    //     A[iSize][iSize], entries are A[row][col]
    // Output:
    //     return value is TRUE if successful, FALSE if pivoting failed
    //     A[iSize][iSize], inverse matrix

    static bool Solve (int iSize, MgcReal** aafA, MgcReal* afB);
    // Input:
    //     A[iSize][iSize] coefficient matrix, entries are A[row][col]
    //     B[iSize] vector, entries are B[row]
    // Output:
    //     return value is TRUE if successful, FALSE if pivoting failed
    //     A[iSize][iSize] is inverse matrix
    //     B[iSize] is solution x to Ax = B

    static bool SolveTri (int iSize, MgcReal* afA, MgcReal* afB, MgcReal* afC,
        MgcReal* afR, MgcReal* afU);
    // Input:
    //     Matrix is tridiagonal.
    //     Lower diagonal A[iSize-1]
    //     Main  diagonal B[iSize]
    //     Upper diagonal C[iSize-1]
    //     Right-hand side R[iSize]
    // Output:
    //     return value is TRUE if successful, FALSE if pivoting failed
    //     U[iSize] is solution

    static bool SolveConstTri (int iSize, MgcReal fA, MgcReal fB, MgcReal fC,
        MgcReal* afR, MgcReal* afU);
    // Input:
    //     Matrix is tridiagonal.
    //     Lower diagonal is constant, A
    //     Main  diagonal is constant, B
    //     Upper diagonal is constant, C
    //     Right-hand side Rr[iSize]
    // Output:
    //     return value is TRUE if successful, FALSE if pivoting failed
    //     U[iSize] is solution

    static bool SolveSymmetric (int iSize, MgcReal** aafA, MgcReal* afB);
    // Input:
    //     A[iSize][iSize] symmetric matrix, entries are A[row][col]
    //     B[iSize] vector, entries are B[row]
    // Output:
    //     return value is TRUE if successful, FALSE if (nearly) singular
    //     decomposition A = L D L^t (diagonal terms of L are all 1)
    //         A[i][i] = entries of diagonal D
    //         A[i][j] for i > j = lower triangular part of L
    //     B[iSize] is solution to x to Ax = B

    static bool SymmetricInverse (int iSize, MgcReal** aafA,
        MgcReal** aafAInv);
    // Input:
    //     A[iSize][iSize], entries are A[row][col]
    // Output:
    //     return value is TRUE if successful, FALSE if algorithm failed
    //     AInv[iSize][iSize], inverse matrix

protected:
    // tolerance for 2x2 and 3x3 system solving
    static MgcReal ms_fTolerance;
};

//----------------------------------------------------------------------------
inline MgcReal& MgcLinearSystem::Tolerance ()
{
    return ms_fTolerance;
}
//----------------------------------------------------------------------------


#endif
